; Collects a set of BOUT++ outputs. Converts into BOUT's XYZT format
;
; August 2008: Updated to read BOUT++ v0.3 output, which is split in X
;              and Y
; ----------------
; UPDATE:
;   - 20160721 by JGCHEN    Email: cjgls@pku.edu.cn
;     o parallelly collecting data, tested on EDISON platform
;       NOTE & WARNING:
;       1. CTRL+ c sometimes can not stop the threads, CTRL+z to stop IDL
;       2. too much threads (>8) may cause freezing up, firstly choosing suitable y_nthreads
;       3. set keyword `noparallel to use single threads, i.e. same as original collect.pro
;       4. Any bug report and suggestions are gladly welcome.
;

FUNCTION in_arr, arr, val

  w = WHERE(arr EQ val, count)
  IF count GT 0 THEN RETURN, 1

  RETURN, 0
END


PRO collect_sub, xind=xind, yind=yind, zind=zind, tind=tind, data=data, dimsize=dimsize, $
      NXPE=NXPE, MXSUB=MXSUB, MYSUB=MYSUB, MXG=MXG, MYG=MYG, indfiles=indfiles, $
      path=path, prefix=prefix, fext=fext, var=var, quiet=quiet
;
; The procedure is called in Function pcollect in this file.
; Please check the main body of pcollect function to get the meaning of options.
;

  IF quiet EQ 0 THEN PRINT, "collect_sub called"

  ndims = N_ELEMENTS(dimsize)
  tsize = tind[1] - tind[0] + 1
  data=[]

  IF ndims EQ 3 THEN BEGIN
  ; could be xyz or txy
      isxyz = 0
      IF dimsize[0] EQ MZ THEN isxyz = 1
  ENDIF

  ndims_indfiles = SIZE(indfiles, /N_DIMENSIONS)
  size_indfiles = SIZE(indfiles, /DIMENSIONS)

  nindf_x = size_indfiles[0]
  IF ndims_indfiles EQ 2 THEN BEGIN
    nindf_y = size_indfiles[1]
  ENDIF ELSE BEGIN
    nindf_y = 1
  ENDELSE
  IF quiet EQ 0 THEN PRINT, "Shape of File array: ", nindf_x, nindf_y

  FOR iy = 0, nindf_y-1 DO BEGIN
    tmp_data=[]
    FOR ix = 0, nindf_x-1 DO BEGIN
      ; get X and Y processor indices
      i = indfiles[ix, iy]
      pe_yind = FIX(i / NXPE)
      pe_xind = i MOD NXPE

      ; output progress of collection in one line
      percent = (iy*nindf_x + ix + 1)/FLOAT(N_ELEMENTS(indfiles))*100.0
      PRINT, FORMAT="(%'\33[1MProcessing: %5.1F \% \33[1A')", percent

      ; get local y range
      ymin = yind[0] - pe_yind*MYSUB + MYG
      ymax = yind[1] - pe_yind*MYSUB + MYG

      xmin = xind[0] - pe_xind*MXSUB
      xmax = xind[1] - pe_xind*MXSUB

      inrange = 1

      IF (ymin GE (MYSUB + MYG)) OR (ymax LT MYG) THEN BEGIN
        inrange = 0 ; out of Y range
      ENDIF

      IF ymin LT MYG THEN ymin = MYG
      IF ymax GE MYSUB + MYG THEN ymax = MYG + MYSUB - 1

      ; check lower x boundary
      IF PE_XIND EQ 0 THEN BEGIN
        ; keeping inner boundary
        IF xmax LT 0 THEN inrange = 0
        IF xmin LT 0 THEN xmin = 0
      ENDIF ELSE BEGIN
        IF xmax LT MXG THEN inrange = 0
        IF xmin LT MXG THEN xmin = MXG
      ENDELSE

      ; check upper boundary
      IF PE_XIND EQ (NXPE - 1) THEN BEGIN
        ; keeping outer boundary
        IF xmin GE (MXSUB+2*MXG) THEN inrange = 0
        IF xmax GE (MXSUB+2*MXG-1) THEN xmax = (MXSUB+2*MXG-1)
      ENDIF ELSE BEGIN
        IF xmin GE (MXSUB+MXG) THEN inrange = 0
        IF xmax GE (MXSUB+MXG) THEN xmax = (MXSUB+MXG-1) ; ndim = 3 seems ERROR???, 2*MXG
      ENDELSE

      ; calculate global indices
      xgmin = xmin + PE_XIND * MXSUB
      xgmax = xmax + PE_XIND * MXSUB

      ygmin = ymin + PE_YIND*MYSUB - MYG
      ygmax = ymax + PE_YIND*MYSUB - MYG

      IF inrange THEN BEGIN
        filename = path+"/"+prefix+"."+STRTRIM(STRING(i),2)+"."+fext
        IF quiet EQ 0 THEN BEGIN
          PRINT, ""
          PRINT, "Reading from "+filename
          PRINT, "PE_X = " + STRTRIM(STRING(pe_xind),2) + $
            " PE_Y = " + STRTRIM(STRING(pe_yind),2)
          PRINT, "Local indices: ["+STRTRIM(STRING(xmin),2)+"-"+STRTRIM(STRING(xmax),2)+"][" + $
            STRTRIM(STRING(ymin),2)+"-"+STRTRIM(STRING(ymax),2)+"]"
          PRINT, "   =>  Global: ["+STRTRIM(STRING(xgmin),2)+"-"+STRTRIM(STRING(xgmax),2)+"][" + $
            STRTRIM(STRING(ygmin),2)+"-"+STRTRIM(STRING(ygmax),2)+"]"
        ENDIF ELSE IF quiet EQ 2 THEN BEGIN
          ; NOTE: can't use WRITEU procedure in new threads.
          ;       the output is random, so it's meaningless
          outinfo = "Reading from "+filename+": ["+STRTRIM(STRING(xmin),2)+"-"+STRTRIM(STRING(xmax),2)+"][" + $
          STRTRIM(STRING(ymin),2)+"-"+STRTRIM(STRING(ymax),2)+"] -> ["+STRTRIM(STRING(xgmin),2)+"-"+STRTRIM(STRING(xgmax),2)+"][" + $
          STRTRIM(STRING(ygmin),2)+"-"+STRTRIM(STRING(ygmax),2)+"]"
          PRINT, FORMAT="(%'\33[1MProcessing: %5.1F\% %s\33[1A')", percent, outinfo
        ENDIF

        ; set the indices
        IF ndims EQ 4 THEN BEGIN
          inds = LONARR(8)
          inds[0] = LONG(tind[0])
          inds[1] = LONG(tind[1])
          inds[2] = LONG(xmin)
          inds[3] = LONG(xmax)
          inds[4] = LONG(ymin)
          inds[5] = LONG(ymax)
          inds[6] = LONG(zind[0])
          inds[7] = LONG(zind[1])
        ENDIF ELSE IF ndims EQ 3 THEN BEGIN
          inds = LONARR(6)
          IF isxyz THEN BEGIN ; [x, y, z]
            inds[0] = LONG(xmin)
            inds[1] = LONG(xmax)
            inds[2] = LONG(ymin)
            inds[3] = LONG(ymax)
            inds[4] = LONG(zind[0])
            inds[5] = LONG(zind[1])
          ENDIF ELSE BEGIN    ; [t, x, y]
            inds[0] = LONG(tind[0])
            inds[1] = LONG(tind[1])
            inds[2] = LONG(xmin)
            inds[3] = LONG(xmax)
            inds[4] = LONG(ymin)
            inds[5] = LONG(ymax)
          ENDELSE
        ENDIF ELSE BEGIN
          inds = LONARR(4)
          inds[0] = LONG(xmin)
          inds[1] = LONG(xmax)
          inds[2] = LONG(ymin)
          inds[3] = LONG(ymax)
        ENDELSE

        ny_loc = ymax - ymin + 1
        nx_loc = xmax - xmin + 1

        handle = file_open(filename)
        d = file_read(handle, var, inds=inds)
        file_close, handle

        size_d = SIZE(d, /DIMENSION)

        ; stack array in x axis: first dimension
        IF ndims EQ 4 THEN BEGIN  ; [t, x, y, z] -> [x, y, z, t]
          rollaxis_d = FLTARR([size_d[1:3], size_d[0]])
          FOR t=0, tsize-1 DO BEGIN
            rollaxis_d[*, *, *, t] = d[t, *, *, *]
          ENDFOR
          tmp_data = [tmp_data, rollaxis_d]
        ENDIF ELSE IF ndims EQ 3 THEN BEGIN
          IF isxyz THEN BEGIN ; [x, y, z]
            tmp_data = [tmp_data, d]
          ENDIF ELSE BEGIN    ; [t, x, y] -> [x, y, t]
            rollaxis_d = FLTARR[size_d[1:*], size_d[0]]
            FOR t=0, tsize-1 DO BEGIN
                tmp_data = [tmp_data, rollaxis_d]
            ENDFOR
          ENDELSE
        ENDIF ELSE tmp_data = [tmp_data, d]
      ENDIF
    ENDFOR
    ; stack array in y axis: second dimension
    data = [[data], [tmp_data]]
    DELVAR, tmp_data
  ENDFOR
  PRINT
END

FUNCTION pcollect, var=var, xind=xind, yind=yind, zind=zind, tind=tind, t_array=t_array, $
    path=path, use=use, quiet=quiet, debug=debug, prefix=prefix, noparallel=noparallel, $
    x_nthreads=x_nthreads, y_nthreads=y_nthreads, thread_limit=thread_limit
; Options:
;   - noparallel    KEYWORD
;       if set this keyword, then collecting data in current single thread
;   - x/y_nthreads  INT
;       number of threads in x/y axis directions
;       WARNING: large threads number (>8) sometimes leads to idl freezing up.
;           Firstly, the y_nthreads is used
;   - thread_limit  INT
;       limit of total threads used. [default: Min([8, ncpu])
;
  IF NOT KEYWORD_SET(thread_limit) THEN thread_limit = 8
  thread_limit = MIN([thread_limit, !CPU.HW_NCPU])
  IF NOT KEYWORD_SET(prefix) THEN prefix="BOUT.dmp"

  IF NOT KEYWORD_SET(debug) THEN BEGIN
    on_error,2  ; If an error occurs, return to caller

    ; catch any error which occurs
    CATCH, errcode

    IF errcode NE 0 THEN BEGIN
      ; error occurred
      PRINT, "Error occurred in collect: "+!ERROR_STATE.MSG

      RETURN, 0
    ENDIF
  ENDIF

  ; first find out how many files there are

  IF NOT KEYWORD_SET(use) THEN use = 0

  IF NOT KEYWORD_SET(path) THEN path = "./data/"
  path = FILE_EXPAND_PATH(path) ; MUST be absolute path in IDL_IDLBridge

  IF KEYWORD_SET(quiet) THEN quiet = 2 ELSE quiet = 1
  IF KEYWORD_SET(debug) THEN quiet = 0

  IF NOT KEYWORD_SET(var) THEN BEGIN
    var = "Ni"
    PRINT, "No variable specified. Trying to read "+var
  ENDIF

  ;; Allow user to specify a single number for the index

  IF KEYWORD_SET(xind) THEN BEGIN
      IF N_ELEMENTS(xind) EQ 1 THEN xind = [xind, xind]
  ENDIF
  IF KEYWORD_SET(yind) THEN BEGIN
      IF N_ELEMENTS(yind) EQ 1 THEN yind = [yind, yind]
  ENDIF
  IF KEYWORD_SET(zind) THEN BEGIN
      IF N_ELEMENTS(zind) EQ 1 THEN zind = [zind, zind]
  ENDIF
  IF KEYWORD_SET(tind) THEN BEGIN
      IF N_ELEMENTS(tind) EQ 1 THEN tind = [tind, tind]
  ENDIF

  SPAWN, "\ls "+path+"/"+prefix+".*", result, exit_status=status

  IF status NE 0 THEN BEGIN
      PRINT, "ERROR: No data found"
      RETURN, 0
  ENDIF

  nfiles = N_ELEMENTS(result)

  IF quiet EQ 0 THEN PRINT, "Number of files found: ", nfiles

  mainfile = result[use] ; File to get most data

  ; Get the file extension
  i = STRPOS(mainfile, '.', /REVERSE_SEARCH)
  fext = STRMID(mainfile, i+1, STRLEN(mainfile)-(i+1))

  ; Select again only the ones with this extension
  SPAWN, "\ls "+path+"/"+prefix+".*."+fext, result, exit_status=status

  nfo = nfiles
  nfiles = N_ELEMENTS(result)

  IF (quiet EQ 0) AND (nfo NE nfiles) THEN PRINT, "Number of files with same type: ", nfiles

  ; get list of variables
  handle = file_open(mainfile)
  var_list = file_list(handle)

  ; Check if the variable requested is available

  IF MAX(STRCMP(var_list, var)) EQ 0 THEN BEGIN
    PRINT, "Variable '"+var+"' not found"
    ; Check if the case is different
    res = STRCMP(var_list, var, /FOLD_CASE)
    w = WHERE(res EQ 1, count)
    IF COUNT GT 0 THEN BEGIN
      var = var_list[w[0]]
      PRINT, "-> Variables are case-sensitive: Using '"+var+"'"
    ENDIF ELSE BEGIN
      file_close, handle
      RETURN, 0
    ENDELSE
  ENDIF

  var_list = STRUPCASE(var_list)
  ; get the version of BOUT++

  IF in_arr(var_list, "BOUT_VERSION") THEN version = file_read(handle, "BOUT_VERSION") ELSE VERSION= 0.2

  IF quiet EQ 0 THEN PRINT, "BOUT++ VERSION: ", version

  ; read parameters
  mxsub = file_read(handle, "MXSUB")
  mysub = file_read(handle, "MYSUB")
  MZ    = file_read(handle, "MZ")
  myg   = file_read(handle, "MYG")
  t_array = file_read(handle, "t_array")
  nt = N_ELEMENTS(t_array)

  IF version GT 0.25 THEN BEGIN
    ; 2D decomposition - need to read NXPE and MXG

    NXPE = file_read(handle, "NXPE")
    MXG  = file_read(handle, "MXG")
    NYPE = file_read(handle, "NYPE")
    NPE  = NXPE*NYPE

    IF NPE LT nfiles THEN BEGIN
        PRINT, "WARNING: More files than expected"
        nfiles = NPE
    ENDIF ELSE IF NPE GT nfiles THEN BEGIN
        PRINT, "WARNING: Some files missing"
    ENDIF

    IF quiet EQ 0 THEN PRINT, "MXG = ", mxg

    NX = NXPE * MXSUB + 2*MXG
  ENDIF ELSE BEGIN
    NX = MXSUB
    MXG = 0
    NXPE = 1
    NYPE = nfiles
  ENDELSE

  NY = MYSUB * NYPE

  IF quiet EQ 0 THEN BEGIN
      PRINT, "Size of the grid: ", NX, NY, MZ
      PRINT, "In each file: ", MXSUB, MYSUB, MZ
  ENDIF

  ; check that var is 4-D

  ndims = file_ndims(handle, var)
  dimsize = file_size(handle, var)
  ; dimension size [z, y, x, t], netcdf reverses the order of indices

  file_close, handle ; Close the result[use] file

  IF NOT KEYWORD_SET(xind) THEN xind = [0, nx-1]
  IF NOT KEYWORD_SET(yind) THEN yind = [0, ny-1]
  IF NOT KEYWORD_SET(zind) THEN zind = [0, mz-2]
  IF NOT KEYWORD_SET(tind) THEN tind = [0, nt-1]

  ; check ranges

  IF xind[0] LT 0 THEN xind[0] = 0
  IF xind[0] GT nx-1 THEN xind[0] = nx-1
  IF xind[1] LT 0 THEN xind[1] = 0
  IF xind[1] GT nx-1 THEN xind[1] = nx-1
  IF xind[0] GT xind[1] THEN BEGIN
    tmp = xind[0]
    xind[0] = xind[1]
    xind[1] = tmp
  ENDIF

  IF yind[0] LT 0 THEN yind[0] = 0
  IF yind[0] GT ny-1 THEN yind[0] = ny-1
  IF yind[1] LT 0 THEN yind[1] = 0
  IF yind[1] GT ny-1 THEN yind[1] = ny-1
  IF yind[0] GT yind[1] THEN BEGIN
    tmp = yind[0]
    yind[0] = yind[1]
    yind[1] = tmp
  ENDIF

  IF zind[0] LT 0 THEN zind[0] = 0
  IF zind[0] GT MZ-2 THEN zind[0] = MZ-2
  IF zind[1] LT 0 THEN zind[1] = 0
  IF zind[1] GT MZ-2 THEN zind[1] = MZ-2
  IF zind[0] GT zind[1] THEN BEGIN
    tmp = zind[0]
    zind[0] = zind[1]
    zind[1] = tmp
  ENDIF

  IF tind[0] LT 0 THEN tind[0] = 0
  IF tind[0] GT nt-1 THEN tind[0] = nt-1
  IF tind[1] LT 0 THEN tind[1] = 0
  IF tind[1] GT nt-1 THEN tind[1] = nt-1
  IF tind[0] GT tind[1] THEN BEGIN
    tmp = tind[0]
    tind[0] = tind[1]
    tind[1] = tmp
  ENDIF

  xsize = xind[1] - xind[0] + 1
  ysize = yind[1] - yind[0] + 1
  zsize = zind[1] - zind[0] + 1
  tsize = tind[1] - tind[0] + 1
  data = []

  IF ndims GT 1 THEN BEGIN ; ndims= 2, 3 or 4
    ; print ranges
    IF quiet EQ 0 THEN BEGIN
      PRINT, "X indices: ", xind[0], " to ", xind[1]
      PRINT, "Y indices: ", yind[0], " to ", yind[1]
      IF ndims EQ 3 THEN BEGIN
        IF isxyz THEN BEGIN
          PRINT, "Z indices: ", zind[0], " to ", zind[1]
        ENDIF ELSE BEGIN
          PRINT, "T indices: ", tind[0], " to ", tind[1]
        ENDELSE
      ENDIF
    ENDIF

    ; determine index range of X,Y-processor
    r_pe_xind = FIX((xind-MXG)/MXSUB) ; xind including x-guard cells
    r_pe_yind = FIX(yind/MYSUB)       ; yind excluding y-guard cells

    ; check up boundary
    r_pe_xind = MIN([[NXPE-1, NXPE-1], [r_pe_xind]], DIMENSION=2)
    r_pe_yind = MIN([[NYPE-1, NYPE-1], [r_pe_yind]], DIMENSION=2)

    ; thread size used in cxx source code
    x_nthreads_src = r_pe_xind[1] - r_pe_xind[0] + 1
    y_nthreads_src = r_pe_yind[1] - r_pe_yind[0] + 1

    ; index range of dump files in which the target data locate
    r_nfile = INTARR(x_nthreads_src, y_nthreads_src)
    FOR iny = r_pe_yind[0], r_pe_yind[1] DO BEGIN
      FOR inx = r_pe_xind[0], r_pe_xind[1] DO BEGIN
        r_nfile[inx - r_pe_xind[0], iny - r_pe_yind[0]] = iny * NXPE + inx
      ENDFOR
    ENDFOR

    ; parallelly collecting
    ; TODO:
    IF KEYWORD_SET(y_nthreads) THEN y_nthreads=y_nthreads ELSE y_nthreads = y_nthreads_src
    IF KEYWORD_SET(x_nthreads) THEN x_nthreads=x_nthreads ELSE x_nthreads = x_nthreads_src
    y_nthreads = MIN([y_nthreads, thread_limit, y_nthreads_src, 8])
    x_nthreads = MIN([x_nthreads, thread_limit/y_nthreads, x_nthreads_src])

    IF (NOT KEYWORD_SET(noparallel)) THEN BEGIN
      oBridge = OBJARR(x_nthreads, y_nthreads)
      PRINT, FORMAT='(%"Creating %I x %I Threads")', x_nthreads, y_nthreads

      ; get path of first 'pcollect.pro'
      script_dir = FILE_WHICH('pcollect.pro')

      IF quiet EQ 0 THEN PRINT, "STARTING MULTI-THREADS ..."
      FOR iny = 0, y_nthreads-1 DO BEGIN ; parallel processing in y axis
        FOR inx = 0, x_nthreads-1 DO BEGIN ; parallel processing in x axis

          ; get dump file range in current thread
          x_down = inx*x_nthreads_src/x_nthreads
          x_up = (inx+1)*x_nthreads_src/x_nthreads-1
          y_down = iny*y_nthreads_src/y_nthreads
          y_up = (iny+1)*y_nthreads_src/y_nthreads-1

          indfiles = r_nfile[x_down:x_up, y_down:y_up]
          IF quiet EQ 0 THEN BEGIN
            PRINT, 'Thread ', iny*x_nthreads+inx
            PRINT, FORMAT='(%"Dump File range: x=[%I, %I], y=[%I, %I]")', x_down, x_up, y_down, y_up
            ; set keyword `OUTPUT to display message from child process
            oBridge[inx, iny] = OBJ_NEW('IDL_IDLBridge', OUTPUT='')
          ENDIF
          IF iny+inx EQ 0 THEN oBridge[inx, iny] = OBJ_NEW('IDL_IDLBridge', OUTPUT='') $
            ELSE oBridge[inx, iny] = OBJ_NEW('IDL_IDLBridge')
          oBridge[inx, iny]->SetVar, 'xind', xind
          oBridge[inx, iny]->SetVar, 'yind', yind
          oBridge[inx, iny]->SetVar, 'zind', zind
          oBridge[inx, iny]->SetVar, 'tind', tind
          oBridge[inx, iny]->SetVar, 'NXPE', NXPE
          oBridge[inx, iny]->SetVar, 'MXSUB', MXSUB
          oBridge[inx, iny]->SetVar, 'MYSUB', MYSUB
          oBridge[inx, iny]->SetVar, 'MXG', MXG
          oBridge[inx, iny]->SetVar, 'MYG', MYG
          oBridge[inx, iny]->SetVar, 'indfiles', indfiles
          oBridge[inx, iny]->SetVar, 'path', path
          oBridge[inx, iny]->SetVar, 'prefix', prefix
          oBridge[inx, iny]->SetVar, 'fext', fext
          oBridge[inx, iny]->SetVar, 'var', var
          oBridge[inx, iny]->SetVar, 'dimsize', dimsize
          oBridge[inx, iny]->SetVar, 'quiet', quiet
          oBridge[inx, iny]->Execute, '!quiet=quiet'  ; disable outputting compilation info
          oBridge[inx, iny]->Execute, ".r " + script_dir
          oBridge[inx, iny]->Execute, "collect_sub, xind=xind, yind=yind, zind=zind, tind=tind, NXPE=NXPE, MXSUB=MXSUB, MYSUB=MYSUB, MXG=MXG, MYG=MYG, indfiles=indfiles, path=path, prefix=prefix, fext=fext, var=var, data=idata, dimsize=dimsize, quiet=quiet", /NOWAIT
        ENDFOR
      ENDFOR

      notdone = 1
      WHILE notdone DO BEGIN
        done=0
        FOR i=0, N_ELEMENTS(oBridge)-1 DO $
          done = done+oBridge[i]->Status()
        IF done EQ 0 THEN notdone=done
      ENDWHILE

      IF KEYWORD_SET(debug) THEN BEGIN
        oBridge[0]->Execute, "help, /FULL, output=hp_doc"
        hp_doc=oBridge[0]->GetVar('hp_doc')
        PRINT, hp_doc
      ENDIF

      ; stack array
      FOR iny=0, y_nthreads -1 DO BEGIN
        tmp_data = []
        FOR inx=0, x_nthreads -1 DO BEGIN
          ; stack array in x axis: first dimension
          tmp_data = [tmp_data, oBridge[inx, iny]->GetVar('idata')]
          oBridge[inx, iny]->Execute, 'DELVAR, idata'
          OBJ_DESTROY, oBridge[inx, iny]
          PRINT, FORMAT="(%'\33[1MThread %I Finished\33[1A')", iny*x_nthreads + inx
        ENDFOR
        ; stack array in y axis: second dimension
        data = [[data], [tmp_data]]
      ENDFOR
      DELVAR, tmp_data
      PRINT
    ENDIF ELSE BEGIN
      ; collecting data in current one thread
      indfiles = r_nfile
      collect_sub, xind=xind, yind=yind, zind=zind, tind=tind, NXPE=NXPE, MXSUB=MXSUB, MYSUB=MYSUB, MXG=MXG, MYG=MYG, indfiles=indfiles, path=path, prefix=prefix, fext=fext, var=var, data=data, dimsize=dimsize, quiet=quiet
    ENDELSE
  ENDIF ELSE BEGIN
    ; 0 or 1-D
    ; Just read the variable from main file
    IF quiet LT 2 THEN PRINT, "Reading from "+mainfile
    handle = file_open(mainfile)
    data = file_read(handle, var)
    file_close, handle
  ENDELSE

  ;IF quiet EQ 2 THEN PRINT, ""

  CATCH, /CANCEL

  RETURN, data
END
